Loudspeaker cabinet design by topology optimization

Using material distribution-based topology optimization, we optimize the bandpass design of a loudspeaker cabinet targeting low frequencies. The objective is to maximize the loudspeaker’s output power for a single frequency as well as a range of frequencies. To model the loudspeaker’s performance, we combine a linear electromechanical transducer model with a computationally efficient hybrid 2D–3D model for sound propagation. The adjoint variable approach computes the gradients of the objective function with respect to the design variables, and the Method of Moving Asymptotes (MMA) solves the topology optimization problem. To manage intermediate values of the material indicator function, a quadratic penalty is added to the objective function, and a non-linear filter is used to obtain a mesh independent design. By carefully selecting the target frequency range, we can guide the optimization algorithm to successfully generate a loudspeaker design with the required bandpass character. To the best of our knowledge, this study constitutes the first successful attempt to design the interior structure of a loudspeaker cabinet using topology optimization.


The 3D model
Consider the cross-section of the loudspeaker setup illustrated in Fig. 2 (left), with dimensions w × h × d , where d is in a direction perpendicular to the plane.The dimensions of the baffle and the output port are w × d and h p × d , respectively.Let Ŵ c b and Ŵ c f denote the back and the front of the transducer's diaphragm, respectively.Similarly, let b and f denote the back and the front chamber, respectively.Let e c be the unit vector in the direc- tion of motion of diaphragm, the negative y-direction in Fig. 2 (left).Moreover, let Ŵ p denote the output port that separates the loudspeaker's interior from its exterior, and Ŵ s all sound-hard walls.

The loudspeaker's interior
We consider time-harmonic linear wave propagation inside the loudspeaker's interior and assume that the acoustic pressure satisfies P(x, t) = ℜ{p(x)e iωt } , where i is the imaginary unit, p the complex pressure amplitude, ω the angular frequency, and t the time.This assumption gives us the Helmholtz equation for p.That is, (1) where c is the speed of sound and = ∇ • ∇ the Laplace operator.
As stated earlier, we assume all the walls to be rigid (sound-hard), which implies the boundary condition The linearized Euler equation gives a relation between the diaphragm velocity u c and the pressure p b inside b , and the pressure p f inside f , where k = ω/c is the wave number, ρ is the air density, and n b and n f are outward directed normals with respect to b and f , respectively.Furthermore, we use mechanical and electric circuit equations to model the electromechanical properties of the transducer.The mechanical model expresses the balance of forces on the speaker diaphragm, where M md is the moving mass of the diaphragm, R ms the mechanical resistance, C ms the mechanical compliance, Bl the force factor, and I the electric current amplitude in the voice coil.
Remark 1 Throughout this article, we omit symbols of measure, such as dŴ or d , in integral expressions, since the type of measure will be clear from the domain of integration.
To complete the transducer model, we use the simple electric circuit illustrated in Fig. 3 (right) given by where the amplifier voltage V is the input to the system, R is the electric resistance, and L is the inductance.

The loudspeaker's exterior
To model the interaction of the output port Ŵ p with the loudspeaker's exterior, we divide Ŵ p into N p = N (3a) www.nature.com/scientificreports/For each frequency under consideration, we compute the impedance matrix Z p 3D column by column by solving N p exterior Helmholtz problems with successive unit excitations on each of the panels.The boundary-element solver in the commercial Pafec VibroAcoustics software is used for this calculation.

The hybrid model
Consider the cross-section of the loudspeaker illustrated in Fig. 3 (left).For computational purposes, we split the interior into two domains at the baffle, an upper box formed by the domain above the baffle, and a lower box, with dimensions w × h l × d , formed by the domain below the baffle.The domain is the full lower box.We assume that all the walls and internal solid structures inside are extruded in the z-direction.Let γ t denote the boundary that separates the two domains, and let γ p denote the output port.
We assume planar symmetry in the acoustic pressure along the z-axis in the air region of .The variations in the direction of the z-axis will be negligible due to the long wavelengths.This assumption allows the use of 2D wave propagation (in the xy-plane), which provides two advantages.First, it is computationally advantageous for use in an optimization loop compared to a full 3D model.Second, the planar symmetry provides construction advantages, as the interior can then be built by placing, say, wooden slabs aligned with the z-axis in the lower box.However, we employ a 3D model for the upper box and the exterior because 3D effects cannot be avoided there.

The loudspeaker's exterior
We assume planar symmetry in the lower box for the acoustic pressure, but we cannot use the same assumption to model the interaction with the exterior.Similarly, as for the 3D model, we pre-compute the acoustic properties by assembling an impedance matrix.We assume the acoustic velocity to be constant on each boundary segment γ p j , j = 1, . . ., N p h and extend each of these into depth-running strips, γ p j × (0, d) on which the acoustic velocity still is assumed to be constant.Now the exterior response of a unit velocity on each strip can be computed using the full 3D boundary-element method.In fact, the response is already available from the matrix Z p 3D by adding all columns corresponding to a particular strip.For each strip on the port, averaging corresponding rows of Z p 3D , we obtain the average pressure response to the applied unit strip velocity.In this way, we obtain an N

The loudpeaker's interior
The upper box We seek to represent the interaction at γ t with the upper box in terms of the acoustic response of the back chamber as well as the coupling to the electromechanical model of the transducer.Hence, in addition to pressures and velocities, also the diaphragm velocity u c (in the negative y-direction) and the voice coil current I are taken into account while computing the response of the upper box.Although it is reasonable to assume planar symmetry and thus carry out 2D calculations in the lower box, we cannot make this assumption above γ t due to the presence of the cylindrically-shaped transducer that creates a local sound field extending throughout the www.nature.com/scientificreports/back chamber.Therefore we make a "hybrid" ansatz and compute in full 3D in the upper box, but average the pressure response in the depth direction as it concerns the effects on the lower box.More precisely, we divide the boundary γ t into line segments γ t j , where j = 1, 2, . . ., N t , extend these into the depth region to obtain strips γ t × (0, d) .Each such strip, in succession, is given a unit velocity amplitude, while the other strips, as well as the diaphragm velocity, are held at zero velocity.In addition, a unit diaphragm velocity is imposed, while the strips are held at a zero velocity.In each of these cases, the acoustic pressure is then computed in the full 3D back chamber.Finally, for each of these excitations, we compute the voice coil current I as well as the pressure response averaged over each of the N t strips.This procedure enables us to obtain an impedance relation of the form where the N t × N t block Z tt represents the acoustic response of the back chamber, the N t × 1 block Z tc , the 1 × N t block Z ct , the 1 × 1 block Z cc represent the interaction with the speaker diaphragm; T , in which u t j is the normal velocity on γ t j , and p γ t j denotes the average pressure on γ t j .To compute this matrix, we set up a full 3D finite element model of the upper box, including a linear electromechanical model of the transducer, in the commercial software Comsol Multiphysics.Matrix Z t is then computed, for each frequency, column by column by exciting the velocity of each strip in succession as well as an excitation of the diaphragm and computing corresponding voice coil current and averaging the pressure response over each strip.

The lower box
Physically, the panels γ p j and γ t j correspond to massless and stiff pistons.The acoustic impedance relation is valid provided that the boundaries γ p and γ t have air on both sides.To ensure this property, we split into three non-intersecting parts, illustrated in Fig. 4, denoted t , d , and p and do not allow material to be placed in t and p .Further, we define material indicator function α such that α = 1 in the air region and α = 0 in the solid region.(In practice, we use α = 1 in the air region and α = ε > 0 in the solid, where ε is a small positive number.)As mentioned earlier, we do not allow material to be placed in t and p .That is, we require α ≡ 1 in both p and t .Wave propagation in the lower box is governed by the following Helmholtz equation for the acoustic pressure,

Remark 2
The model in Eq. ( 9), where the material indicator function α controls the presence and absence of material inside , was introduced by Wadbro and Berggren 15 for acoustic horn optimization and has been used in many contributions since then.In material distribution topology optimization, redefining the material indicator function by replacing α = 0 with α = ε > 0 is a standard strategy 6 to obtain a unique solution to the governing equation.A relevant question is how much this approximation will affect the solution.Using this approximation in the Helmholtz type Eq. ( 9), Kasolis et al. 24 performed an analysis which revealed that the error is linear in ε to the solution of an exactly modeled scatterer with Neumann conditions on the sound-hard walls.In the computations reported below, we set ε = 10 −3 .( 8) www.nature.com/scientificreports/

Finite element discretization
The domain is discretized using a uniform mesh of square elements.We use standard continuous, piecewise bi-quadratic basis functions for the complex pressure, denoted by ϕ 1 , ϕ 2 , . . ., ϕ N � , where N is the number of degrees of freedom in the finite element approximation.In addition, we approximate α by an element-wise constant material indicator function α h .By putting together a finite element discretized version of the governing Eq. ( 9), impedance relations ( 8) and (7), and circuit Eq. ( 5), we arrive at the following equation system where the entry p j in vector p represents the complex pressure amplitude on the jth degree of freedom, and the N × N matrix A ll , the N × N p matrix A lp , the N × N t matrix A lt , the N p × N matrix A pl , and the N t × N matrix A tl have entries respectively, and finally a Ic = Bl/(ρc) and a II = R + iωL.

Design definition by filtering
To enable gradient-based optimization, we relax α h to allow intermediate values that neither represent solid material nor air.We employ a combination of a penalization and an appropriate filtering approach to enforce extreme values of α h and ensure size control.In this section, we detail our definition and handling of the design variables using a non-linear filtering method.The morphological dilate and erode operators are approximated by these non-linear filters.More precisely, to use gradient-based optimization, we approximate the morphological operators with the harmonic-mean based filters suggested by Svanberg and Svärd 25 , and we use the fW-mean based filtering framework of Wadbro and Hägg 26 .
For the optimization, d is our design domain.We define the set of admissible design variables as where d is a N d × 1 vector that defines the material distribution inside d before filtering.
For the filtering, we use the extended domain f = d ∪ s ∪ a illustrated in Fig. 4 (right).Here, s is a domain occupied by solid, and a = t ∪ p is a domain occupied by air.The uniform mesh of square elements is extended to f using elements of the same size as those in .The N elements of f are sorted so that those in d come first, s second, and a last.We assemble a weight matrix W r = D −1 r G r , where G r is a neighborhood indicator matrix with entries and D r = diag(G r 1 N ) , where 1 N = (1, 1, . . ., 1) T ∈ R N .In expression (13), x i and x j are the centroids of elements i and j in f , respectively, and x i − x j is the distance between x i and x j .In addition, we define the functions where β > 0 is a parameter.We denote the inverse functions of f E β and f D β by f −1 E β and f −1 D β , respectively.Now, we define the discrete harmonic erode and dilate operators that act on an N × 1 vector η with entries 0 ≤ η N ≤ 1 as www.nature.com/scientificreports/respectively.Here, T .(To keep the discussion simple, we refer to the harmonic-mean based approximate morphological operators as harmonic followed by the name of the morphological operator that they approximate.)Parameter β governs the properties of the filter, which approaches a linear blurring filter in the limit β → +∞ .The non-linearity of the filter increases as β decreases.We thus will refer to β as the non-linearity parameter.In the limit β → 0 , the action of the filters tends to that of the corresponding morphological operators.Using the harmonic dilate and erode operators in a series, we define the harmonic close operator on N d × 1 vector d as where I N d is the identity matrix of size N d × N d .More precisely, in expression ( 16), first we expand the design to the extended domain f , then we apply the filter on f , and finally we extract the filtered entries of d .Finally, we define (cf.Remark 2) the vector F (d) , which holds the element values of α h in d , as The optimization problem The radiated power of the loudspeaker through the output port is the integral of the Poynting vector over the port.In the discretized case, the radiated power becomes where the overline denotes complex conjugate, x is the solution to Eq. ( 10), x * is the Hermitian transpose of x , and, using the same blocking as in Eq. ( 10), A suitable way to evaluate the performance of this type of loudspeaker is to assume the loudspeaker to be placed on an infinite floor in an anechoic chamber and measure the sound pressure level (SPL) in dB at 1 m in front of the output port, for a given input voltage.Since the SPL is proportional to the logarithm of the radiated power, we base our design objective on this quantity.More precisely, to optimize for a set of frequencies f 1 , f 2 , . . ., f m , we define the objective function where P f i is the output power evaluated according expression (18) in the case where x solves governing Eq. ( 10) for angular frequency ω = 2πf i and physical design α h with element values F (d) , as defined in expression (17), in d .To solve the optimization problem, we use the method of moving asymptotes (MMA) by Svanberg 27

Sensitivity Analysis
The MMA algorithm requires the gradients of the objective function with respect to design variables.We employ the adjoint variable method because of its efficiency; it computes the full gradients at the cost of (at most) one extra finite element analysis.Here, we present the sensitivity analysis for the radiated power given by expression ( 18). ( (  www.nature.com/scientificreports/Perturbing α h and using that B is real and symmetric, we obtain a first order perturbation of the radiated power given by Similarly, the first order perturbation of Eq. ( 10) is Pre-multiplying Eq. ( 24) with an arbitrary vector z * , we obtain Now we select z as the solution to the so-called adjoint equation which yields that x * B = z * A , which substituted into expression (23) gives where last equality follows from expression (25).
Recall that α h is element-wise constant.Letting n be the nth element in d -so a n is the element value of α h in n -expression (27) implies that where (28) provides sensitivities with respect to the physical design a , which in turn depends on design variables d , which are the ones updated by the MMA algorithm, through a compound function involving the non-linear filter, as detailed in expression (17).Using the chain rule, as presented in detail by Hägg and Wadbro 29 together with expression (28), we obtain the gradient of P with respect to d.

Numerical experiments and results
Consider a loudspeaker box (recall Fig. 3 (left)) with dimensions w × h × d = 80 cm × 100 cm × 60 cm .A baffle of dimensions w × d = 80 cm × 60 cm , containing an 18-inch transducer, separates the upper and the lower box.The lower box has dimensions w × h l × d = 80 cm × 70 cm × 60 cm .An output port of dimensions h p × d = 30 cm × 60 cm is located at the lower left of the loudspeaker.The electromechanical parameters of the transducer, given in Table 1, are typical for a commercial 18-inch driver.
Two types of numerical experiments are performed.First, the loudspeaker is optimized for single frequencies ranging from 20 to 100 Hz .Second, it is optimized for selected frequency bands.The hybrid model is implemented in MATLAB, and all experiments are performed using a finite element discretization with 280 × 320 square elements, which yields a side length of each element of 2.5 mm .With this resolution, we solve for 363,502 unknowns in Eq. (10).
We start the optimization with d = (1, 1, . . ., 1) T ∈ R N d as initial guess, that is, with no solid material inside the design domain d .We employ a continuation technique, in which the penalty and non-linearity parameters are gradually increased and decreased, respectively.Initially, the material distribution problem is solved with a very small value of the penalty parameter ( ζ = 10 −4 ).The MMA algorithm computes the residual norm of the Karush-Kuhn-Tucker (commonly referred to as the KKT or the first-order optimality) conditions.We increase the penalty parameter and use the previous solution as the initial guess when the KKT residual norm is less than 10 −2 times the initial KKT residual norm.The intermediate values of d become expensive as the penalty increases with each step.As a result, when the penalty parameter is sufficiently large, all entries of d are essentially 0 or 1 www.nature.com/scientificreports/ at optimum.For the non-linear filter, we start with a higher value of the non-linearity parameter ( β = 10 ) and decrease the non-linearity parameter whenever we increase the penalty parameter.Throughout the optimization, the relation ζ = 10 −3 β −1 holds.
To evaluate the performance of the loudspeaker, we use SPL 1 m .Assuming that the loudspeaker is placed on an infinite sound-hard floor, the quantity SPL 1 m is defined as where p 1m is the pressure at the floor 1 m in front of the output port, and p o = 20 µPa is the reference pressure amplitude for computing SPL.
Furthermore, we consider the layout with an empty lower box as our reference loudspeaker, which is also the initial guess for the optimization.For all the results, we compare the SPL 1 m of all the optimized designs (single as well as multi-frequency optimization) with the SPL 1 m of the reference loudspeaker.To compute the SPL 1 m , we use V = 1 V as the input voltage to the amplifier in all cases.

Single-frequency optimization
The resulting designs and corresponding frequency responses for single-frequency optimizations using target frequencies ranging from 20 to 80 Hz are shown in Fig. 5.The optimized loudspeakers' frequency response (solid blue line) is compared to the reference loudspeaker's frequency response (dashed black line).The reference curve has a peak at 60 Hz .We note that the frequency response of the optimized design yields a peak at the frequency subject to optimization.Here, the optimization algorithm tunes the resonance frequency of the system to the frequency subject to optimization.Figure 5a,b show that the peak SPL 1 m of the designs for the lower target frequencies (20 and 40 Hz ) is lower than the peak SPL 1 m of the reference loudspeaker.However, the peak SPL 1 m of the designs for the higher target frequency (80 Hz ) is close to that of the reference curve (Fig. 5d).As we already have a peak at 60 Hz , the optimization algorithm for this target frequency does not put solid material (Fig. 5c) inside the lower box, and the frequency response overlaps the reference curve.This suggests that the initial design (empty lower box) is a local minimum to the optimization problem for 60 Hz .The improvement in the efficiency at the target frequency is between 0 dB (for 60 Hz ) and +10 dB (for 80 Hz ).At the lowest frequency in the range, 20 Hz , the peak of frequency response is at 86.5 dB .This corresponds to an increase in efficiency of approximately +11.5 dB at 20 Hz compared to the empty box.Here, the single-frequency optimization is only used for preliminary investigation.The main objective of this study is to perform multi-frequency optimization of the loudspeaker enclosure presented in the following section.

Multi-frequency optimization
Here, we optimize the loudspeaker over single-octave bands, one-and-a-half octave bands, and a double-octave band.The optimization frequencies are logarithmically spaced within the design frequency band.That is, f i = F start 2 i/12 , for i = 0, 1, 2, . . ., n − 1 , with F start being the lowest frequency in the range, and n = 13 , n = 19 , or n = 26 if the frequency band is a single octave wide, one-and-a-half octave wide, or a double octave wide, respectively.Here, the optimization algorithm adds very little solid material inside the lower box; that is, the optimization shows that the empty (initial) design is close to a local minimum for the targeted frequency band.However, for the last two bands, which are, 45-90 Hz and 50-100 Hz , the frequency responses in Fig. 6c,d show that there are improvements in the output power over the full frequency band subject to optimization.
The frequency bands for the one-and-a-half octave band optimization are 30-85 Hz , 32.5-92 Hz , 35-99 Hz , and 37.5-106 Hz .Fig. 6e-h show the results.The frequency response overlaps the reference curve in Fig. 6e, and there is no improvement in performance.Again, the optimized design is very similar to the initial design, with very little solid material within the design domain.This indicates that the initial design is close to a local minimum.To avoid this minimum, we chose a slightly higher starting frequency and considered the frequency band 32.5-92 Hz .The frequency response in Fig. 6f demonstrates an improved bandpass design with a larger, more distinct pass band compared to the 30-85 Hz band.The results for 35-99 Hz and 37.5-106 Hz , shown in Fig. 6g,h, are qualitatively similar.
For the double-octave band optimization, the frequency band is 30-120 Hz , and the results are shown in Fig. 7.In comparison to the one-and-a-half octave band, the frequency response shows a wider bandpass design and slightly better performance for low frequencies than the reference curve.Comparing Fig. 6c with 6g and 6d with 6h, the single octave band results look superficially better than the one-and-a-half octave results.This is due to the higher SPL peak achieved inside the target frequency range for the latter, which compensates for the lower levels outside the target frequency band.The one-and-a-half octave band objective function is indeed worse for the single-octave-band design, compared over the corresponding one-and-a-half octave band.Above the design frequency band, the SPL curve can behave unpredictably, which is in analogy with for instance filter design or high order polynomial interpolation.

Single-frequency optimization
As previously stated, by placing solid material into the design domain, the optimization algorithm tunes the system's resonance frequency to the frequency under consideration.It should be noted that tuning the resonance frequency to the optimization target frequency is accomplished in two ways, depending on whether the target frequency is below or above the resonance frequency for the reference design.In addition to the results presented in Fig. 5, we performed additional experiments targeting other frequencies.For all our experiments, optimizing for lower frequencies yields a reduction in the port size while optimizing for higher frequencies results in a reduction in the lower box volume.
The lower box acts as a Helmholtz resonator 30 because the wavelength is large compared to the dimensions of the lower box.The standard lumped model for the resonance frequency of the Helmholtz resonator is (2π) −1 c A p /(L p V l ) , where A p and L p are the area and effective acoustic length of the output port, respectively, and V l is the volume of the lower box.According to Fig. 5, there is more material inside the lower box for 20 Hz as compared to 40 Hz, and reducing the volume of the lower box increases the resonance frequency, and vice versa.It seems that a slightly increased opening could be compensated by a reduction in volume.It should be noted that the single-frequency case tends to invoke extreme designs when the target frequency is pushed toward the physically lower limit.The designs of Fig. 5a,b contain narrow channels which would likely induce nonlinear effects for appreciable sound levels.If, for some reason, a single frequency subwoofer would be of practical interest, a design constraint precluding too narrow channels would be advisable.

Multi-frequency optimization
To achieve a bandpass design, a dividing wall appears in the lower box for the single octave bands 45-90 Hz and 50-100 Hz , as well as for the one-and-a-half-octave bands 32.5-92 Hz , 35-99 Hz , 37.5-106 Hz , and the double octave-band 30-120 Hz .The design can be seen as a cascade of two Helmholtz resonators.A lumped model analysis suggests that this loudspeaker box corresponds to a sixth-order acoustic filter with a high-frequency roll-off of 36 dB per octave.This is the limit behavior occurring for frequencies well above the design frequency band.The even faster roll-off seen for the presented designs is governed by the quality factor of the Helmholtz resonators.The lumped filter cascade model assumes that the filter elements somehow retain their acoustical identity when connected.This assumption can be far from accurate due to the near field interaction between the resonators in the subwoofer and can be interpreted as one culprit for the failure of the lumped model.

Validation of the hybrid model
Recall that for computational efficiency, we use a hybrid model, where we assume planar symmetry in the acoustic pressure along the z-axis in the lower box.It is reasonable to ask how accurate this assumption is, particularly when solid material is placed inside the lower box.To assess this issue, we consider the design optimized for the one-octave design (45-90 Hz ), depicted in Fig. 6c.For this design, we use COMSOL Multiphysics to imple- ment a full 3D model, as detailed by Bokhari et al. 23 , with second-order tetrahedral elements with a maximum side length of 0.05 m.The 3D model uses a body-conforming mesh and sound-hard boundary conditions at the interfaces between air and solid material in the lower box.Figure 8 shows the loudspeaker's frequency response computed using the hybrid and the 3D model.This result validates the fidelity of the hybrid model.

Figure 1 .
Figure 1.Cutaway drawing of the loudspeaker with an empty front chamber.

, �p� Ŵ p 2 ,
panels, Ŵ p j , where j = 1, 2, . . ., N p , as illustrated in Fig.2(right).For each panel, we let u p j denote the complex normal velocity on Ŵ p j and p Ŵ p j the average pressure on Ŵ p j .We assemble an N p × N p matrix Z p 3D that relates the normal velocities u p = u p 1 , u p 2 , . . ., u p N p T to the pressures p p = �p� Ŵ p 1 . . ., �p� Ŵ p N p T through the impedance relation (2) ∂p ∂n = 0, on Ŵ s .

Figure 2 .
Figure 2. The 3D model.Left: Cross-section of the loudspeaker with an empty front chamber.Right: Front view of the loudspeaker where Ŵ p is divided into N p = N p h × N p d panels.
the normal velocities u p on each strip γ p j × (0, d) to the average pressures p p on all strips through the impedance relation where u p = u

Figure 3 .
Figure 3. Left: The hybrid model, cross-section of loudspeaker with an empty front chamber.Right: The electric circuit model of the transducer.

Figure 5 .
Figure 5.Single-frequency optimization: The figures above the graphs show the material distribution inside the domain after optimization.The boundary at the top indicates the 18-inch transducer, and the boundary at the lower left is the output port.In the graphs, the dashed black line is the frequency response of the reference loudspeaker (empty lower box), and the solid blue line is the frequency response of the optimized system for the corresponding frequency.

Figure 6 .
Figure 6.Multi-frequency optimization for single and one-and-a-half octave bands: The figures above the graphs show the material distribution inside the domain after optimization.The boundary at the top indicates the 18-inch transducer, and the boundary at the lower left is the output port.In the graphs, the dashed black line is the frequency response of the reference loudspeaker (empty lower box), and the solid blue line is the frequency response of the optimized system for the selected frequency range.

Figure 7 .
Figure 7. Multi-frequency optimization for double-octave band: Left: Figure shows the material distribution inside the domain after optimization.The boundary at the top indicates the 18-inch transducer, and the boundary at the lower left is the output port.Right: In the graph, the dashed black line is the frequency response of the reference loudspeaker (empty lower box), and the solid blue line is the frequency response of the optimized system for the selected frequency range.

Figure 8 .
Figure 8.The frequency response for the optimized design in Fig. 6c evaluated using the hybrid model (solid blue line) and the 3D model (red diamonds).
. Since the MMA expects a minimization problem, we write the optimization problem as By relaxing the material indicator function optimization, we obtain physical designs with intermediate values.As previously stated, this is undesirable; thus we deal with intermediate values using penalization and filtering.To suppress the intermediate values in d , we add a quadratic penalty term 28 to the objective function, which results in the problem where ζ is a positive penalty parameter.

Table 1 .
Electromechanical properties of the 18 inch transducer.